function [ SSt,AAt,PPt] = KalmanSmoother( YY,TT,RR,QQ,DD,ZZ,HH,ZZ_t,msel,varargin )

% /**********************************************************
% ** Kalman filter
% **
% ** {loglh,obserror,obsvar,AT,PT} = kalmanf(YY,TT,RR,QQ,DD,ZZ,HH,VV,ZZ_t,A0,P0,tcond);
% **
% **  y(t) = DD + ZZ*x(t) + ZZ_t*t+ u(t)
% **  u(t) ~ N(0,HH)
% **
% **  x(t) = TT*x(t-1) + RR*e(t)
% **  e(t) ~ N(0,QQ)
% **
% **  cov(e(t),u(t)) = VV
% **********************************************************/
loglhzero = -1E8;

nobs   = size(YY,1);
ny     = size(YY,2);

obserror  = zeros(nobs,ny);
obsvar    = zeros(nobs,ny);
loglh     = 0;

if nargin<=9
     scale_initial_var=0.1;
     At=zeros(size(TT,1),1);
     Pt=dlyap(squeeze(TT(:,:,1)),squeeze(RR(:,:,1))*squeeze(QQ(:,:,1))*squeeze(RR(:,:,1))'); %??????????????????????????????

else
    At = varargin{1};
    Pt = varargin{2};
end

%% Filter: Going forward
t = 1;
AAt=zeros(size(TT,1),nobs);
PPhat=zeros(size(TT,1),size(TT,1));
PPt=zeros(size(TT,1),size(TT,1));
while t <= nobs
    
    % 		/*trend and constant */
    DD = DD + ZZ_t;
    
    %Adjust the matrix ZZ in case of missing observations    
    ZZm=ZZ(~isnan(YY(t,:)),:);
    DDm=DD(~isnan(YY(t,:)));
    YYm=YY(t,~isnan(YY(t,:)));
    At1 = At;
    Pt1 = Pt;
    
    % 		/** Forecasting **/
    TTm=squeeze(TT(:,:,msel(t))); RRm=squeeze(RR(:,:,msel(t))); QQm=squeeze(QQ(:,:,msel(t)));
    alphahat = TTm*At1;
    Phat     = TTm*Pt1*TTm' + RRm*QQm*RRm';
    Phat     = .5*(Phat+Phat');

    
    
    yhat     = ZZm*alphahat + DDm;
    nut      = YYm - yhat';
%     YYm
%     yhat'
%     keyboard
    
    
    Ft       = ZZm*Phat*ZZm' + HH(~isnan(YY(t,:)),~isnan(YY(t,:)));% + ZZ*RR*VV + (ZZ*RR*VV)';
    Ft       = .5*(Ft+Ft');
    
    
%     if t > tcond;
        loglh    =  loglh -0.5*size(YYm,2)*log(2*pi) - 0.5*log(det(Ft)) - 0.5*nut*pinv(Ft)*nut';
%     end;

    At = alphahat + Phat*ZZm'*pinv(Ft)*nut';
    Pt = Phat - (Phat*ZZm')*pinv(Ft)*(Phat*ZZm')';
    
    %store the estimated state
    AAt(:,t)=At;
    PPhat(:,:,t)=Phat;
    PPt(:,:,t)=Pt;
    
    t = t+1;
    
end
SSt=[];

%% Smoother: Going backward
SSt=zeros(size(TT,1),nobs);
SSt(:,nobs)=AAt(:,nobs); 
PPT(:,:,nobs)=PPt(:,:,nobs);
for t=nobs-1:-1:1
    TTm=squeeze(TT(:,:,msel(t+1))); RRm=squeeze(RR(:,:,msel(t+1))); QQm=squeeze(QQ(:,:,msel(t+1)));
  [PPhatInv]=PseudoInverse(PPhat(:,:,t+1));
    SSt(:,t)=AAt(:,t)+PPt(:,:,t)*TTm'*PPhatInv*(SSt(:,t+1)-TTm*AAt(:,t));
    PPT(:,:,t)=PPt(:,:,t)+(PPt(:,:,t)*TTm'*PPhatInv)*(PPT(:,:,t+1)-PPhat(:,:,t+1))*(PPt(:,:,t)*TTm'*PPhatInv)';
end
end


